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Abstract 

This paper reviews the essential physics of gravitational instability in a Robertson- Walker back- 
ground spacetime. Three approaches are presented in a pedagogical manner, based on (1) the Eu- 
lerian fluid equations, (2) the Lagrangian description of trajectories, and (3) the Lagrangian fluid 
equations. Linear and nonlinear limits are discussed for each case. Shear and tides are shown to 
play a key role in nonlinear gravitational instability. 

The Lagrangian fluid approach is used to show that several widely held beliefs about gravita- 
tional instability are false. The following collapse theorem is proven: for a given initial density 
fluctuation and growth rate, the spherical tophat perturbation collapses more slowly than any 
other conflguration. We also show that density maxima are not the flrst points to collapse and 
that underdense regions may collapse if their initial shear is sufficiently high. The Lagrangian fluid 
approach leads to an almost closed set of local evolution equations for individual mass elements. 
The magnetic part of the Weyl tensor, which may be present even in the nonrelativistic (Newto- 
nian) limit, may prevent a purely local description. However, neglecting the magnetic part of the 
Weyl tensor, we obtain predictions for high-redshift collapse that are in good agreement with a 
high-resolution cold dark matter N-body simulation. 



1 Introduction 

Gravitational instability theory provides useful relations between the large-scale mass (or galaxy) 
density field p{x) and velocity field v(^x). These relations allow us, in principle, to achieve the 
following goals: 

1. Reconstruct three-dimensional fields from distances and redshifts {Hor, cz) ([3]) 

2. Recover initial fluctuation fields ([26], [15]) 

3. Test theories of initial fields ([18], [11], [30], [28]) 

4. Test gravitational instability paradigm ([10]) 



5. Determine cosmological parameters ([19], [10]) 



This program constitutes a growing subfield of cosmological research. 

Achieving these goals requires use of theoretical tools of gravitational instability. This paper 
summarizes some of the methods and recent results, with particular emphasis on the Lagrangian fluid 
method. 



2 Eulerian fluid equations 

The best-known dynamical description of mass is given by the Eulerian fluid equations. We express 
them here using the comoving position x and conformal time r, which are related to the proper 
position r and time t measured in a locally inertial comoving frame by 

g = r/a{t) , dr = dt/a{t) , (1) 

where a(f) is the cosmic expansion scale factor. Because a > 0, r is a monotonic function of t, so that 
we may write a = ^(t). (E.g., for an Einstein-de Sitter universe, a oa t^l^ implies r a t^l^ and a a r^.) 
We also refer to the density fluctuation field 6{x,t) and "peculiar" velocity field v(^x,t), which are 
related to the proper density p and velocity df/dt as follows: 

The reader will recognize H as the Hubble parameter and should be familiar with its relation to the 
mean density p and the density parameter Q = 8TrGp/{3H^). 

The density and velocity fields are assumed to evolve consistently with conservation of mass and 
momentum. For a nonrelativistic perfect fluid on scales much smaller than the Hubble distance c/H , 
this implies the fluid equations ([5]): 

Continuity: — + V • [(1 + ^)'i;] = , (3) 

OT 



Euler: -^^ + (v -V^v = -- v - V (j) 

OT a 



1^ 

-Vp 

P 



(4) 
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Poisson: = ^T^Ga^Sp = - QoH^a-^S . (5) 

These are identical in form with the noncosmological fluid equations aside from the Hubble drag term 
in the Euler equation (N.B. d/a = aH is not the Hubble parameter!) and the fact that the mean 
density p does not contribute to the gravitational potential (j). (In the literature (j) is often called 
"peculiar," but this is carrying historical precedent too far. There is nothing peculiar about either (j) 
or V.) We assume that nongravitational (e.g., electromagnetic) forces are unimportant on large scales. 

The fluid equations are valid for individual components (e.g., coUisionless dark matter or baryons) 
as well as for the combined distribution of mass as long as the stress tensor is isotropic (i.e., given by 
the pressure p). Before trajectories intersect (or, it is believed, after the fluid variables are spatially 
averaged so that 8^ ^1)) the fluid equations with p = provide a good approximation to the dynamics 
of coUisionless dark matter. For baryons, particle collisions are sufficiently rapid so that the stress is 
isotropic in the fluid frame and the pressure term must be included. Terms depending on pressure are 
kept in brackets in the following so that the reader can readily see how pressure modifies the evolution. 

A key consequence of equation (4) is that gravity alone cannot generate vorticity uj = V X v. 
Taking the curl of equation (4) we obtain 



duj -> d 

— = S/ x(vxuj) UJ 

OT a 



-^{Vp) X (Vp) 
P 



(6) 



The expression in brackets is called the baroclinic term. It vanishes for a pressureless or an isentropic 
fluid (i.e., one with constant specific entropy) or, more generally, for a barotropic fluid with p = p{p). 



/ 



dp 
P 



(8) 



The baroclinic term is very important in meteorology (it generates cyclones, tornados, and hurricanes!) 
but not cosmology. From equation (6) it follows that irrotational = 0) flow remains irrotational in 
the absence of baroclinic torques (Kelvin's circulation theorem). 

Equation (6) implies that any primeval vorticity decays as w a unless it is so large that 
turbulent stresses amplify vorticity through the v X uj term. The latter possibility seems to be ruled 
out by the isotropy of the cosmic microwave background radiation. In this case, the velocity field 
prior to the epoch of galaxy formation (and today, on large scales) is expected to be a potential flow 
(i.e., irrotational). Potential flow is the essential assumption underlying the POTENT analysis ([3]). 
If w = 0, we may write the velocity field in terms of a potential $: 

v=-V^, ^{x,t) = ^{0,t)- r v-dl. (7) 

Jo 

Any path may be used to evaluate the velocity potential $, in particular, a radial path for which 
$(r, e, 4>) - $(0) = - /o (cz' - Hor') dr'. In the POTENT procedure ([3], [9], [4]), we first interpolate 
the redshifts and distance of galaxies in the local universe to define a smooth radial velocity field. The 
radial integral gives the velocity potential, whose transverse derivatives give the components of the 
smoothed velocity field not obtainable directly from redshifts and distances. 

For a potential flow, the Euler equation may be reduced to an evolution equation for the velocity 
potential: 

5$ 1 , ,o d 
Bernoulli: V$r = ^ + (l> 

The Bernoulli, continuity, and Poisson equations are the basis for several approximate nonlinear 
solution methods ([25], [15], [20], [23], [1]). 

The Eulerian fluid equations are commonly linearized to give the evolution at early times or on 
large scales such that <C 1. Assuming a potential flow, the velocity field is described fully by its 

— > 

divergence, 6 = V ■ v. Linearizing equation (3) and the divergence of (4), and assuming isentropic 
variations with = {dp/ dp)s'^ p = P (where S is the specific entropy and Cg is the sound speed), 
we get 

1^ + ^ = 0, ^ = -^e-A'KG-paH-clvH . (9) 

OT OT a 

These two equations may be combined to give a second-order (in time) linear differential equation for 
6. The spatial dependence is solved using plane waves with comoving wavenumber k. It is easy to see 
that gravitational (Jeans) instability results if k^c^ < ^'KGpa? . For = (or on any scales that are 
strongly Jeans unstable), the general solution for the density and velocity is a linear combination of 
growing and decaying modes D±{t): 

6 = 6+{x)D+{t) + 6-{x)D_{t) , 9 = -6+{x)D+{t)- 6-{x)D4t) . (10) 

At late times the growing mode dominates. Its logarithmic derivative with respect to the expansion 
factor depends on the background cosmology chiefly through Q, and is written din D^/ din a = /(fi) ~ 
([27], [22]). 

Sometimes it is stated that 9 <x 6 solely as a consequence of the continuity equation. This is 
false. If the linear growing mode dominates (and if the pressure, vorticity, and nongravitational 
forces are negligible), then 9 = —aHf{Q,)6. However, if the decaying mode is present or any of 
the other conditions is violated, then 9 is no longer proportional to 6, despite mass conservation. 
Because the growing mode quickly overtakes the decaying mode, unless a perturbation is created with 
unnaturally small ^_|_/^_, the growing mode should dominate by the present time so that independent 
measurements of 6 and 9 may be used to estimate fijl). (This is true even if the perturbations were 
created by nongravitational processes, provided that gravity subsequently drives D-/D^ 0.) Dekel 
et al. ([10]) have used a quasi-nonlinear generalization of this idea to place bounds on Q. 

In the linear regime (i.e., while eqs. 9 are valid), the velocity and density of individual Fourier 
components evolve independently. When the perturbations become large different harmonics are 
coupled. The evolution is no longer local either in Fourier or real space. However, the evolution 
becomes somewhat easier to follow if we abandon the Eulerian description for a Lagrangian one. 



3 Lagrangian description of trajectories 



In a Lagrangian description one follows individual mass elements, in contrast with the Eulerian practice 
of tracking the values of the fluid variables at fixed spatial coordinates. Different mass elements (or 
particles) are labeled by a fixed Lagrangian coordinate q, so that the trajectories are x = x(^q,T). If 
there are no forces except gravity, the motion is governed by Newton's laws in comoving coordinates: 

d'^x ( d'^x\ d dx -> , , , 

— > 

This equation is to be solved for each mass element after we relate — to the trajectories. How? 

— > 

In general, — must be computed by solving the Poisson equation with the mass distribution 
given by the instantaneous positions of all the mass elements. There are, however, several circum- 
stances in which this computation simplifies. 

First, if the mass distribution has sufficient symmetry, the Poisson equation may be solved using 
Gauss' theorem. For example, if the mass distribution is spherical about a point ^ = 0, the solution is 

- V</. = — ^ [M(r) - M(r)] , (12) 

where M(r) = (47r/3) /o(ar)'^, r = \x\, and eV = x/r. Before trajectories intersect, M is fixed for a 
given mass element and is therefore a Lagrangian coordinate. Equations (11) and (12) together give an 
equation of motion for r(T, M) for fixed M. The equation can be solved to give exact radial Keplerian 
trajectories. Identical results are obtained whether one uses Newton's laws in a noncosmological 
background or full general relativity ([27], §§19, 87). 

Gauss' theorem may also be applied in the case of planar symmetry (and, of course, cylindrical 
symmetry). Rather than carrying out this derivation, we consider an alternative approach due to 
Zel'dovich, which is exact for plane-parallel perturbations of cold dust (coUisionless matter with no 
thermal velocities) and provides a good approximation for small perturbations of arbitrary geometry. 

First we choose q (without loss of generality) so that x(^q,0) = q. Then, in general, x(^q,T) = 
q + i/j(^q,T), where i/j(^q,T) is called the displacement field. Mass conservation implies 



S{x,t) _ ^ 



<^(^,0) 



streams 



d 



dq 



(13) 



where the sum is taken over all the q present at a given x. In the linear regime there is only one stream 
and the perturbations to the Jacobian determinant are small. Then, following Zel'dovich ([32]), we 
expand the Jacobian to first order in di/j''/dx^. Assuming ^(^, 0) = 0, we get the following relation 
between 6 and 6 = {4TrGa^p)~^ V • ~ — Vg* - which yields (to first order in 

-A'KGa^p'if . (14) 

Substituting equation (14) into equation (11), one obtains a linear second-order ordinary differential 
equation in time for i/j(^q,T), whose general solution is 

i;{q,T) ^ D+{T)i;+{q) + D_{T)iP_{q) . (15) 

The same functions D±{t) appear here as in equation (10). 

Zel'dovich proposed extending equation (15) into the nonlinear regime. Given the trajectories, it 
is straightforward to obtain the velocity and density fields as a function of the Lagrangian position q, 

— * 

although more work is required to express them in terms of the Eulerian position x = q — ip. It is well 
known that the Zel'dovich approximation works very well until trajectories intersect and mass elements 
collapse to infinite density; it is superior to the linear Eulerian treatment in which mass elements 
never collapse. However, the Zel'dovich approximation breaks down after trajectories intersect. It is 



equivalent to a one timestep N-body algorithm in that particles are pushed in a fixed direction by 
an amount proportional to the initial acceleration, with no allowing for trajectories to reverse after 
crossing potential minima. There are various ways to cure this problem, including adding viscosity 
(converting the Euler equation to Burgers' equation), pre-filtering the density field, and using higher- 
order perturbation theory (equivalent to taking more than one timestep). We refer the interested 
reader to the review article of Shandarin & Zel'dovich ([29]) as well as to several contributions in this 
volume (and [7], [16], [21], [14]) for further discussion of the Zel'dovich approximation and extensions. 



4 Lagrangian fluid equations 

The Lagrangian approach may be applied to give equations of motion for fluid variables in addition 
to trajectories. We will derive the Lagrangian fluid equations for a pressureless gas beginning from 
the Eulerian equations. The same results follow for cold dust by considering trajectories of adjacent 
freely-falling mass elements. 

The procedure we follow is to rewrite the fluid equations using the Lagrangian (or convective) time 
derivative 

Applied to any Eulerian field /(^, r), df/dr gives the time derivative following the fluid since the fluid 
velocity is i; = dx/dr. 

Replacing the Eulerian time derivatives by Lagrangian ones, the zero-pressure fluid equations in 
comoving coordinates become ([6], but note that the present definition of uj differs by a factor of 2) 

Continuity: — + {1 + 6)6 = , (17) 

ar 

Raychaudhuri: ^ + -9+- 9^+ a'^a.j - - u)^ = -A-KGa^pS , (18) 

ar a 3 2 

Vorticity: ^ + - + - 61 - C7\- w^' = , (19) 

ar a 3 

Shear: ^'^^'■^^\^ ^"-^ o-^fc^-^j- + ^ ^i^j - ^ ^ij (^(^''^(^kl + = -Eij ,(20) 

where we have defined the shear and tide tensors: 

(Note that repeated indices are to be summed over and that we are implicitly assuming the use of 
Cartesian comoving coordinates.) The shear and tide tensors are symmetric and traceless. 

It is remarkable that, aside from Eij, equations (17)-(20) provide a closed set of local evolution 
equations for 6, 9, iJ, and aij, with no spatial derivatives aside from those implicit in the tide tensor. 
Thus, in the absence of gravity, each mass element evolves independently of all the others, at least 
until it crosses other elements (recall that we have neglected pressure). Note that in the linear regime, 
equations (17) and (18) reduce to (9) (neglecting pressure) because d/dr ~ d/dr. However, the 
Lagrangian equations have an advantage in retaining a local form in the nonlinear regime. 

The locality of the fluid equations inspires us to go further to try to obtain an evolution equation 
for Eij. In the Newtonian framework this is unnatural (though not illegal!): gravity is given instan- 
taneously by solution of the Poisson equation and not by some time evolution equation. However, 
in general relativity it is natural to treat the tide tensor on the same footing as the fluid variables. 
Ellis ([12], [13]) derives the appropriate equations using the Bianchi identities and Einstein equations. 
Written in comoving coordinates, the evolution equation for Eij is 

^ + 1 E,, + 9E,, - 3a\^E^)k + ct^'^Em + \ e^^ - Vfc e^^ ^E^-^i = -A7rGa^p{l + 6) a,, . (22) 



Parentheses indicate symmetrization: a^^-Ej^^. = | [a^^Ejk + cr^jEik). The fully antisymmetric Levi- 
Civita tensor is eijk (with £123 = +!)• Equation (22) can be derived (with some difficulty) as an exact 
equation in the Newtonian limit ([17]). The difficult part is the term involving the symmetric traceless 
tensor Hij, called, in general relativity, the magnetic part of the Weyl tensor. Ellis says that Hij has 
no Newtonian analogue, but that is not completely correct, although its Newtonian interpretation 
at present remains unclear. On the other hand, the tide tensor Eij, known in general relativity as 
the electric part of the Weyl tensor, is easily understood in the Newtonian framework. The Weyl 
tensor is the traceless part of the Riemann tensor (the tensor responsible for Newtonian tidal forces). 
The electromagnetic terminology is used because Eij and Hij obey equations similar to the Maxwell 

— * — > — * — * 

equations. For example, equation (22) is analogous to the Ampere law dE/dt + V X -B = AttJ. 

Aside from the term involving Hij, equation (22) is purely local. Thus, if Hij = 0, we have 
obtained a closed set of local Lagrangian equations for the nonlinear evolution of cold dust. This fact 
was first noted by Barnes & Rowlingson ([2]) and applied in cosmology by Matarrese, Pantano, & 
Saez ([24]). Bertschinger & Jain ([6]) fully explored the consequences of this assumption for nonlinear 
evolution of cold dust in an Einstein-de Sitter universe. All of these authors assumed that if w = 
and aij oc Eij initially, then Hij = at least until trajectories intersect. The physical interpretation 
is that local Newtonian evolution is causal and exact in the weak-field limit if mass motions do not 
generate gravitomagnetism or gravitational radiation. 

However, recent analytical and numerical results suggest that Hij in equation (22) does not vanish 
identically in the Newtonian limit. The algebraic complexity of the equations makes it difficult to 
analyze the behavior of the Weyl tensor in the Newtonian limit. Nevertheless, neglecting Hij may 
provide a good approximation in many circumstances, and we will follow the consequences of this 
assumption in Section 7 below. 



5 Four propositions 

Consider the evolution of irrotational, pressureless matter under gravity. Would you agree with the 
following four propositions? 

Proposition A: For a given 6 and 6, a spherical tophat perturbation is the configuration that collapses 
most rapidly. 



Proposition B 



Proposition C 



Proposition D 



Initial density maxima are the sites where collapse first occurs. 

Underdense regions do not collapse before colliding with other streams. 

The final stage of collapse is generically one-dimensional, leading to Zel'dovich pan- 



cakes. 

Now consider four alternative propositions. 

Proposition 1: For a given 6 and 6, a spherical tophat perturbation is the configuration that collapses 
most slowly. 

Proposition 2: Initial density maxima are not the sites where collapse first occurs. 

Proposition 3: Underdense regions do collapse before colliding with other streams, 2/ the initial shear 
is not too small. 

Proposition 4: The final stage of collapse is generically iiyo-dimensional, leading to strongly prolate 
filaments. 

Propositions 1-3 are true while A-C are false. Moreover, if the magnetic part of the Weyl tensor 
is negligible. Proposition 4 is true also. In the following we will see show how Propositions 1-4 follow 
from the Lagrangian fluid equations. 



6 Collapse theorem and corollaries 



Combining the Lagrangian continuity and Raychaudhuri equations we obtain the following exact 
equation for cold dust: 

= ^-^^ + {l + S)(a''a,,-lu;' + 4'KGa'ps) . (23) 
a 6 1 + \ I J 

This equation makes no assumptions about w or B.ij\ it is independent of the Weyl tensor. We can 
use it to investigate gravitational collapse oo) of mass elements before their trajectories intersect 

others. 

If w = 0, equation (23) shows that collapse is accelerated by nonzero shear {a'^^aij is non-negative). 
The spherical tophat perturbation, with 6{x,Ti) = Si for r < R and for r > R, has uniform radial 
flow for r < R, with i; a so that aij = 0. Therefore, for a given 6 and 6, the spherical tophat 
perturbation collapses more slowly than any other configuration (Proposition 1, now a Theorem). The 
physical explanation is that shear increases the rate of growth of the convergence of fluid streamlines. 

Corollary 1: If ^ > 0, an irrotational growing-mode perturbation collapses in finite time ([6]). 

Corollary 2: The collapse time depends on all three initial eigenvalues of d^cj)/ dx^dx^ (equivalently, 
for growing mode perturbations, / dx^dx^), not just on the trace (i.e., 6). Therefore, local maxima 
of 6 are not, in general, local minima of the collapse time (Proposition 2). 

Corollary 3: If aij 7^ 0, a perturbation may collapse even if ^ < initially (Proposition 3). This 
follows from a continuity argument: If ^ = initially, a zero-shear perturbation barely avoids collapse 
in finite time. Nonzero shear induces a positive 6, speeding up collapse, which then occurs in a finite 
time. 

Propositions 1-3 also follow under the Zel'dovich approximation, where the density depends on 
all three eigenvalues of the strain tensor dx^/dq^. However, the Zel'dovich approximation used to 
evaluate this tensor is only an approximation, whereas equation (23) is exact. 

7 Results of local evolution 

Proposition 4 is difficult to analyze because the geometry of collapse depends on the shear tensor, 
whose evolution depends on the tide. In the Zel'dovich approximation, collapse is purely kinematical: 
because the displacements are proportional to the initial accelerations, collapse occurs first in the 
direction corresponding to the smallest eigenvalue of the initial strain tensor. (Not all points will 
collapse because in some places all three eigenvalues oi dip^ j dq^ will be positive.) However, this is not 
exact because the Zel'dovich approximation neglects the gravitational feedback of the collapsing flow 
on itself. 

Bertschinger & Jain ([6], see also [8]) have analyzed equation (22) and shown that the —3a^^-Ej^j^ 
shear-tide coupling term drives the evolution toward prolate configurations of the shear and tides. If 
Hij = 0, the shear tensor generically becomes strongly prolate (with two negative and one positive 
eigenvalue) as collapse is approached. Note that this does not imply that the instantaneous shape of 
a collapsing object is prolate, because dx^/dq^ is related to the time integral of the velocity gradient 
tensor. Nevertheless, the result is surprising because it contradicts the Zel'dovich pancake paradigm. 

There is a plausible physical explanation for prolate collapse. The gravitational binding energy is 
larger for a linear mass distribution (with ^ a Inr ^ —00 for a line mass) than for a planar mass 
distribution (^ constant as r ^ for a sheet of mass). Although the gravitational energy is still 
larger for a one- dimensional (spherical) collapse, nonzero shear and angular momentum conservation 
prevent spherical collapse in general. 

Assuming Hij = 0, Bertschinger & Jain also showed that for a growing-mode Gaussian random 
initial density field in an Einstein-de Sitter universe, 56% of initially underdense cold dust mass 
elements collapse. Only 22% of the mass can avoid collapse before intersecting other mass elements. 
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Figure 1: Top: particles in one Lagrangian slice of a high-resolution cold dark matter N-body simula- 
tion at times specified by the specified linear amplitude as- The Eulerian x and y comoving positions 
are shown for particles in a single Lagranian layer = constant. Note how filamentary the structures 
are; the extent in the third dimension is small. Bottom: smoothed contours of constant initial 6 and 
collapse time, the latter computed assuming that the magnetic part of the Weyl tensor vanishes. 



Recently, S. White has investigated numerically the collapse of a uniform ellipsoidal overdense 
perturbation embedded in an Einstein-de Sitter universe. His numerical results agree well with ap- 
proximate analytical theory neglecting the tidal feedback with the surroundings ([31]) and disagree 
quantitatively with the predictions of the Bertschinger-Jain theory assuming Hij = 0. Assuming there 
is no subtle error in the simulation, one is forced to conclude Hij ^ 0. This is difficult to check 
analytically, and the status of Hij remains unclear, although it is known that the Hij term in equation 
(22) vanishes in cases of high symmetry. However, in the general case it seems appropriate to regard 
the neglect of Hij as an approximation whose range of validity requires further investigation. 

The Hij = predictions are tested in Figure 1 against a high-resolution (20 kpc comoving force 
softening distance) N-body simulation of the standard cold dark matter model with 288'^ particles in a 
100 Mpc cube. The first structure to collapse in the simulation was found at high redshift (cg = 0.05 
corresponds to z = 20 if the model is normalized to the COBE quadrupole). The structure in Figure 
1 lies nearly in the x-y plane, indicating that the collapse is strongly filamentary, as predicted. At 
a later time (cg = 0.075) the strongest filament breaks into two clumps, which subsequently merge. 
Contours of constant linear 6 are shown in the bottom left panel; if collapse occurs when 6 = 1.686 
as predicted by the spherical tophat model, then the 6 contours correspond to collapses occurring 
at 1 + Zc = (12,14,16,18,20). (These are underestimates because the density field was smoothed 
slightly to produce the contours.) The lower right panel shows the results for 1 + Zc (with contour 
levels 12,14,16,18,20,22,24) computed by solving equations (17)-(22) assuming Hij = 0. (The ag 
labels apply only to the top two panels and not the bottom ones.) Note that both sets of contours are 
plotted in Lagrangian space; they purport to show (under the approximations of the spherical tophat 
model or Hij = 0) the Lagrangian volumes that should have collapsed by a given redshift. The point 
of initial collapse is rather close to the density maximum in Lagrangian space (in this case Proposition 
B is nearly correct), but both are shifted relative to the point of Eulerian collapse by long wavelength 
displacements. One can see that the spherical tophat model predicts a later collapse (supporting 
Proposition 1, the Collapse Theorem), with less material having collapsed by a given redshift (because 
of its neglect shear), and it does not suggest the formation of filaments. 

This test is not very precise, but it shows that the local description makes some predictions that 
are close to what happens in this N-body simulation. Given the high resolution of the simulation, the 
strong small-scale power of the cold dark matter spectrum, and the strong nonlinearity of the density 
distributions. Figure 1 must be counted a real success of the Lagrangian fluid description. 

8 Conclusions 

Gravitational instability theory still has surprises waiting to be uncovered. Eulerian theory is well 
developed but rather stale, while studies based on Lagrangian trajectories remain fertile. The recent 
applications of Lagrangian fluid dynamics (not SPH!) to cosmology have begun to attack fully nonlinear 
problems that previously could be addressed only in cases of high symmetry. The completeness of 
this approach remains unclear, however, as long as the Newtonian (and relativistic!) behavior of Hij 
(the magnetic part of the Weyl tensor) is not understood. It seems remarkable that basic questions 
about Newton's laws still exist, but this fact only serves to demonstrate the richness of nonlinear 
gravitational clustering. 

Lagrangian methods are not a panacea. Even if the local description of the fluid variables prior 
to trajectory-crossing were to prove correct, the description is incomplete without knowledge of the 
Eulerian location of the mass elements. This fact suggests that it may prove fruitful to combine 
Lagrangian fluid dynamics with integration of the trajectories to supplement, or even replace ([24]), 
standard N-body techniques. 
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